Disclaimer 1: The present material is an adaptation of the original R script prepared by Dr. Matteo Fontana for the a.y. 2020/2021 Nonparametric statistics course, and later Prof. Andrea Cappozzo (a.y 2022-2023) While I acknowledge Matteo and Andrea for the (great) work done I hereby assume responsibility for any error that may be present in this document.
Disclaimer 2: I will start from the assumption that you are all intermediate R users, if this is not the case please let me know, and we shall find a solution together.
library(MASS)
library(rgl)
library(DepthProc)
library(hexbin)
library(aplpack)
library(robustbase)
library(MDBED) # plot exponential bivariate distribution
Let us start by simulating \(500\) bivariate datapoints whose distribution is marginally exponential
set.seed(2781991) # reproducibility
n=500
df_bivariate_exp = cbind(rexp(n), rexp(n))
head(df_bivariate_exp)
## [,1] [,2]
## [1,] 0.3004277 0.35894610
## [2,] 2.4430839 0.22392999
## [3,] 1.8444128 0.15064919
## [4,] 0.2383688 1.12290176
## [5,] 1.9613656 0.07966151
## [6,] 5.2480885 0.01762207
and visualize their scatterplot
plot(df_bivariate_exp[,1],df_bivariate_exp[,2], xlab="exp 1", ylab="exp 2")
We can further employ a hexagonal binning plot to visualize the data density
bin=hexbin(df_bivariate_exp[,1],df_bivariate_exp[,2], xbins=10, xlab="exp 1", ylab="exp 2")
plot(bin, main="Hexagonal Binning")
Now on depths: there are many possible packages hosted on CRAN to
work with depths in a multivariate setting. We will hereafter use
DepthProc, which is a good general purpose package and in
addition it has good plotting capabilities. Even though not all depth
measures introduced in class are directly available (simplicial, Oja,
convex hull peeling depth), it includes the two depth measures we will
be using throughout this section, namely Tukey and
Mahalanobis depths. Let us look at the help file for
the depth function.
Calculating depth for a given dataset is immediately accomplished by typing
set.seed(20)
tukey_depth=depth(u=df_bivariate_exp,method='Tukey')
It may be useful (spoiler alert, it will become apparent why when we study onparametric forecasting) to calculate the depth of a point relative to a sample. You can do it by:
depth(u = c(0.5, 0.5), X = df_bivariate_exp, method = 'Tukey')
## Depth method: Tukey
## [1] 0.254
Compute the median (deepest point) with depthMedian
function
depthMedian(df_bivariate_exp,depth_params = list(method='Tukey'))
## [1] 0.8364350 0.9297927
Or, having already computed the Tukey depth for the entire sample:
df_bivariate_exp[which.max(tukey_depth),]
## [1] 0.8364350 0.9297927
If you have the luck of dealing with a bivariate dataset,
you can easily visualize the depth surface in a very convenient way.
DepthProc offers you two possible methods:
DepthProc::depthContour(df_bivariate_exp,depth_params = list(method='Tukey'))
or
depthPersp(df_bivariate_exp,depth_params = list(method='Tukey'))
For additional special effects:
depthPersp(df_bivariate_exp,depth_params = list(method='Tukey'),plot_method = 'rgl')
Since the Tukey depth is obtained by counting points in the space of all hyperplanes(surfaces) in \(\mathbb{R}^d\), when \(d>2\), the computational burden increases. Libraries usually recur to heuristic methods which lead to poor estimates (we will see this in the functional data depth measures lab)
a = depthPersp(df_bivariate_exp,depth_params = list(method='Tukey'))
The very same analysis can be carried out by considering a different depth measures:
Simplicial (or Liu) depth:
mrfDepth::sdepth(df_bivariate_exp)$depthZ
## [1] 0.112827245 0.043805925 0.052095903 0.099119154 0.025286042 0.006000000
## [7] 0.038261101 0.146898085 0.008824782 0.189377599 0.141395176 0.015819446
## [13] 0.015610691 0.230319965 0.006000000 0.144052394 0.235806746 0.154630176
## [19] 0.178300939 0.144672139 0.041322742 0.052061617 0.006516310 0.006000000
## [25] 0.064843373 0.035115774 0.151183041 0.032581259 0.168724533 0.075996040
## [31] 0.076487240 0.014556776 0.017281648 0.215428978 0.077749620 0.087816307
## [37] 0.228942318 0.149266485 0.231516141 0.228164812 0.021211387 0.091886037
## [43] 0.174609122 0.031402033 0.009396093 0.137307482 0.047747785 0.013995847
## [49] 0.064876500 0.136080933 0.141716493 0.202258493 0.038173793 0.162751189
## [55] 0.167804525 0.156815511 0.022217688 0.011608470 0.022862399 0.038947340
## [61] 0.125680228 0.071467224 0.048179395 0.150906536 0.024388391 0.009300529
## [67] 0.028471207 0.006373760 0.065631649 0.010840766 0.010340343 0.086454065
## [73] 0.087399039 0.018473574 0.084852645 0.029433952 0.014314750 0.051906077
## [79] 0.040759398 0.142870174 0.087026970 0.112785668 0.200215032 0.025810271
## [85] 0.013132723 0.220927832 0.053925152 0.062973948 0.092953039 0.019528310
## [91] 0.175128136 0.020190453 0.083810416 0.247232054 0.027455151 0.093649757
## [97] 0.131388898 0.037729483 0.113098534 0.058061472 0.032074945 0.007998358
## [103] 0.034047179 0.017662892 0.006189294 0.026123234 0.029044692 0.127519231
## [109] 0.088771422 0.043270927 0.121903132 0.227215733 0.025504455 0.015007026
## [115] 0.217532318 0.185444914 0.010363619 0.036335273 0.068984234 0.029277398
## [121] 0.020245551 0.194301760 0.043730256 0.068852162 0.036250042 0.153333366
## [127] 0.174292633 0.021759905 0.177407490 0.141731077 0.144373663 0.157634208
## [133] 0.105141367 0.112218992 0.079795446 0.046371007 0.036009561 0.154733032
## [139] 0.128767994 0.030082816 0.075373735 0.046179540 0.009522611 0.065790859
## [145] 0.215487553 0.131210904 0.108150276 0.014290364 0.125446314 0.142174614
## [151] 0.091941908 0.237619287 0.234031871 0.149322162 0.210289398 0.134839317
## [157] 0.009364995 0.042450829 0.112996692 0.097948910 0.076924548 0.159482531
## [163] 0.009424729 0.020996547 0.032876355 0.014878432 0.094244730 0.043231137
## [169] 0.018420842 0.009918681 0.175285945 0.006648043 0.006000000 0.037928242
## [175] 0.192059010 0.096792621 0.028786295 0.006823671 0.082096192 0.196624671
## [181] 0.018690586 0.203099645 0.062382983 0.130567448 0.203177584 0.089904097
## [187] 0.063008620 0.013599826 0.047570466 0.025388609 0.132268151 0.223749620
## [193] 0.010616124 0.181931574 0.146335997 0.188846705 0.026875921 0.112860565
## [199] 0.008519931 0.166297656 0.075174445 0.010163315 0.050812130 0.238632784
## [205] 0.008371683 0.032559094 0.089776227 0.007795060 0.188754569 0.006000000
## [211] 0.219730401 0.012509211 0.193274887 0.078380810 0.132731777 0.091073859
## [217] 0.164242895 0.009730835 0.226219089 0.067238622 0.090459859 0.015671681
## [223] 0.019475143 0.007183620 0.036713910 0.148235217 0.052316682 0.022229664
## [229] 0.149962383 0.045797040 0.166040611 0.006071710 0.008082623 0.100052491
## [235] 0.087995557 0.234032064 0.047979960 0.006095420 0.011830746 0.066110776
## [241] 0.131828138 0.013506483 0.038720139 0.156290460 0.007411981 0.014224111
## [247] 0.056601685 0.228979549 0.065634884 0.211382089 0.182183065 0.080841104
## [253] 0.029272811 0.032994712 0.043482435 0.035512181 0.024525919 0.018416302
## [259] 0.212896057 0.035479827 0.043720356 0.246222276 0.056309921 0.157564720
## [265] 0.246196200 0.234541951 0.150499408 0.063904387 0.225730449 0.023263249
## [271] 0.010770505 0.010604293 0.121982568 0.032744380 0.068786151 0.112709564
## [277] 0.126749209 0.131819350 0.038298235 0.029107227 0.172471642 0.144846367
## [283] 0.011919115 0.035989425 0.059580172 0.018198518 0.016025545 0.244591786
## [289] 0.123340512 0.227583504 0.014720477 0.090419586 0.124363909 0.033742231
## [295] 0.169045851 0.009196658 0.173329985 0.009098969 0.172314074 0.038714972
## [301] 0.082954149 0.243680325 0.222993747 0.006000000 0.019619818 0.101167298
## [307] 0.206802328 0.011461574 0.011817563 0.022908274 0.130585412 0.045820895
## [313] 0.006000000 0.042339040 0.072620228 0.171963590 0.012823913 0.217507159
## [319] 0.131091436 0.119595722 0.063843011 0.028823430 0.015158268 0.211571094
## [325] 0.006634764 0.245146631 0.222121061 0.114795326 0.023690176 0.168820050
## [331] 0.175471425 0.190391723 0.030776831 0.088089142 0.246577058 0.119702731
## [337] 0.017343844 0.202930970 0.176901369 0.132935027 0.010365985 0.009306855
## [343] 0.010799334 0.212724051 0.162729797 0.026289446 0.023675930 0.038878383
## [349] 0.013983485 0.050880991 0.170724727 0.031582104 0.065803607 0.078078036
## [355] 0.112293116 0.067284304 0.171121424 0.006000000 0.006000000 0.063771447
## [361] 0.228675906 0.193270686 0.198568993 0.132339329 0.108548326 0.110324746
## [367] 0.057557959 0.113764010 0.124622691 0.009608760 0.016955212 0.016568511
## [373] 0.008966173 0.171838569 0.009098100 0.225907091 0.017107516 0.020021006
## [379] 0.075964169 0.021267933 0.083915349 0.093347080 0.217454330 0.204529348
## [385] 0.102622353 0.008511433 0.184174421 0.006213970 0.204944202 0.201729193
## [391] 0.096672381 0.006000000 0.084770312 0.097965425 0.042973224 0.044008209
## [397] 0.073059082 0.026855977 0.220022696 0.248786537 0.071409856 0.197489775
## [403] 0.109334718 0.070016756 0.058649830 0.047266002 0.032400753 0.081144892
## [409] 0.021029674 0.204261197 0.078311418 0.007849337 0.006000000 0.096135162
## [415] 0.011256682 0.175142816 0.101600116 0.083376874 0.065253640 0.224461212
## [421] 0.006000000 0.011608567 0.038358017 0.032649299 0.077289132 0.057689403
## [427] 0.013298549 0.140689234 0.015727069 0.006402685 0.174051911 0.241760147
## [433] 0.006000000 0.205702924 0.148022358 0.014803197 0.015973489 0.048148152
## [439] 0.025040153 0.088875389 0.073833933 0.026484342 0.015265905 0.204336432
## [445] 0.037277302 0.204826279 0.103561822 0.041502861 0.241818094 0.225796267
## [451] 0.097488664 0.006047903 0.020918222 0.036242075 0.094477533 0.047743197
## [457] 0.149372915 0.058438419 0.074095082 0.065148466 0.008410218 0.031138035
## [463] 0.194908854 0.203296038 0.032386315 0.129884444 0.173873772 0.060589661
## [469] 0.029493010 0.111200135 0.139690658 0.184330299 0.057656808 0.007309800
## [475] 0.228035203 0.168880315 0.055882898 0.055419755 0.020802666 0.144387281
## [481] 0.028750320 0.080912669 0.243760485 0.107519666 0.085550281 0.122074704
## [487] 0.138308907 0.071798006 0.118166212 0.032619214 0.037605235 0.106827148
## [493] 0.087675206 0.122298766 0.025643190 0.068792766 0.092982688 0.007552551
## [499] 0.102314364 0.009563851
and now Mahalanobis depth
maha_depth <- depth(df_bivariate_exp,method='Mahalanobis')
This can be easily hard coded, recall the definition you have seen in class:
\[{M}_{h}D({F;X} ^ {n}) = \frac{ 1 }{ 1 + {{(x - \mu_F)} ^ {T}}{{\Sigma} ^ {-1}}(x - \mu_F) }\]
sample_mean <- colMeans(df_bivariate_exp)
sample_S <- cov(df_bivariate_exp)
maha_depth_manual <- 1/(1+mahalanobis(x = df_bivariate_exp,center = sample_mean,cov = sample_S))
And check that the obtained result is equal to the one obtained via
the depth function ( beware of comparing
floats)
all(abs(maha_depth-maha_depth_manual)<1e-15) # food for thought: sqrt(2) ^ 2 == 2?
## [1] TRUE
all.equal(as.vector(maha_depth)+1, as.vector(maha_depth_manual))
## [1] "Mean relative difference: 0.6672736"
depthMedian(df_bivariate_exp,depth_params = list(method='Mahalanobis'))
## [1] 0.9211997 1.0228290
df_bivariate_exp[which.max(maha_depth_manual),]
## [1] 0.9211997 1.0228290
And again the graphical outputs:
depthContour(df_bivariate_exp,depth_params = list(method='Mahalanobis'))
depthPersp(df_bivariate_exp,depth_params = list(method='Mahalanobis'))
depthPersp(df_bivariate_exp,depth_params = list(method='Mahalanobis'),plot_method = 'rgl')
Please note that anything that comes out of depthContour
or depthPersp is NOT a density: we are not doing
density estimation here, we are performing data exploration by means of
a nonparametric procedure. Let us plot the probability density function
(NOT the depth function) of the bivariate exponential distribution we
drew a sample from:
PDF_3dPlot(rho = 0, Betax = 1, Betay = 1, title="PDF (NOT depth function) of the bivariate exponential distribution", zlabel="")
depthPersp(cbind(rexp(1e4), rexp(1e4)),depth_params = list(method='Tukey'), plot_title="Tukey depth for the exp distribution")
We notice that the tukey depth concentrates around the median, where as the exponential bivariate function increases as x,y go to 0.
Try out the routines seen so far on something even more exotic:
set.seed(1992)
df_bivariate_cauchy = cbind(rcauchy(n,location=0,scale=.001), rcauchy(n,location = 0,scale=.001))
head(df_bivariate_cauchy)
## [,1] [,2]
## [1,] -0.0024184800 0.0008798105
## [2,] -0.0006282804 0.0009932694
## [3,] -0.0006975300 0.0010347342
## [4,] -0.0038080204 0.0141419885
## [5,] 0.0061098712 0.0014929490
## [6,] -0.0012696290 0.0008865863
The Cauchy distribution is a distribution that has neither the first moment (the mean) nor the second moment (the variance). This means that the CLT does not apply.
In the previous Section we have computed depth measures for a given
dataset and we have seen how to effectively plot them. This resulted in
a convenient way to nonparametrically explore a (bivariate) dataset.
Nevertheless, one of the main aims for a statistician to employ depth
measures is to perform multivariate outlier detection. To appreciate
this, let us simulate \(100\) data
points, of which \(95\%\) comes from a
multivariate normal with mean vector mu_good and covariance
matrix sigma_common and the other 5% from another
multivariate normal, which we assume is our outlier generator
process.
mu_good = c(0,0)
mu_outliers = c(7,7)
sigma_common = matrix(c(1,.7,.7,1), ncol = 2)
frac = .05
n=100
# sample points
n_good=ceiling(n*(1-frac))
n_outliers=n-n_good
df_contaminated_normals = data.frame(rbind(
mvrnorm(n_good, mu_good, sigma_common),
mvrnorm(n_outliers, mu_outliers, sigma_common)
))
Let us visualize the true nature of our dataset
class <- c(rep(1,n_good),rep(2,n_outliers))
plot(df_contaminated_normals,xlab="Norm 1", ylab="Norm 2",col=class)
We can clearly see those red points up there… how can we flag them in an automated fashion? The depth contour plot surely helps
depthContour(
df_contaminated_normals,
depth_params = list(method = 'Tukey'),
points = TRUE,
colors = colorRampPalette(c('white', 'navy')),
levels = 10,
pdmedian = F,
graph_params = list(cex=.01, pch=1),
pmean = F
)
But, as we have seen in class, a very handy graphical tool for
spotting multivariate outliers is the bagplot. The bagplot
function from the aplpack package can be used to display
bagplots for bivariate data
bagplot(df_contaminated_normals, factor = 3, show.whiskers = T)
Remember:
bagplot function, we see
that we have a lot of room for customization:aplpack::bagplot(df_contaminated_normals,show.whiskers = F,main="Bagplot")
aplpack::bagplot(df_contaminated_normals,show.loophull = F,main="Sunburst plot")
In addition, if we save the output of the bagplot to an object, we can automatically extract the outliers
bagplot_cont_normals <- bagplot(df_contaminated_normals)
outlying_obs <- bagplot_cont_normals$pxy.outlier
Once the outlying units have been identified, one can discard them from the original data and keep working on the clean subset only. There are several ways to do this.
A more “sql-oriented” approach (this is just a code alternative):
df_clean_1 <-
dplyr::anti_join(
x = df_contaminated_normals,
y = outlying_obs,
by = c("X1" = "x", "X2" = "y"),
copy = TRUE
)
A more “object-oriented programming” approach:
ind_outlying_obs <- which(apply(df_contaminated_normals,1,function(x) all(x %in% outlying_obs)))
df_clean_2 <- df_contaminated_normals[-ind_outlying_obs,]
all.equal(df_clean_1,df_clean_2)
## [1] TRUE
Data for the Hertzsprung-Russell Diagram of the Star Cluster CYG OB1, which contains \(47\) stars in the direction of Cygnus, from C.Doom.
data(starsCYG, package = "robustbase")
names(starsCYG)
## [1] "log.Te" "log.light"
plot(starsCYG, main="Star Cluster CYG OB1")
We can see two groups of points: the majority which tends to follow a steep band, the so called Main Sequence, and four stars in the upper-left corner. In astronomy the \(43\) stars are said to lie on the Main sequence and the four remaining stars are the red giants, namely points with indexes 11, 20, 30 and 34. In details, the red giants are very bright, but they emit light with a very low color-temperature (and thus their surface temperature is still fairly low). We can easily isolate them thanks to the procedure seen so far. Let us look at the perspective plot of the depth surface
depthContour(as.matrix(starsCYG), depth_params = list(method='Tukey'), points=TRUE)
As you can see, the sample mean is biased due to the presence of the 4 red giants, whereas the Tukey Median is not. Let us compute it:
depthMedian(starsCYG)
## log.Te log.light
## 4.38 5.02
As before, we can use the bagplot to visualize and flag the outlying data points.
bagplot_starsCYG <- with(starsCYG,aplpack::bagplot(log.Te,log.light))
red_giants <- bagplot_starsCYG$pxy.outlier
ind_outlying_obs <- which(apply(starsCYG,1,function(x) all(x %in% red_giants)))
ind_outlying_obs
## [1] 11 20 30 34
So far, we have only dealt with bivariate data. Even though the computational complexity escalates quickly when it comes to compute depth measures in \(R ^d\), with \(d>2\), we can still appreciate their usefulness in moderate dimension. Let us generate a trivariate dataset with contamination.
mu_good = rep(0,3)
mu_outliers = c(12,12,3)
sigma_common = diag(3)*2
frac = .1
n=300
# sample points
n_good=ceiling(n*(1-frac))
n_outliers=n-n_good
df_3 = data.frame(rbind(
mvrnorm(n_good, mu_good, sigma_common),
mvrnorm(n_outliers, mu_outliers, sigma_common)
))
class <- c(rep(1,n_good),rep(2,n_outliers))
pairs(df_3, col=class)
To visualize the outliers in this context we retort to a bagplot matrix:
bagplot_matrix <- aplpack::bagplot.pairs(df_3)
You can notice that outlier detection becomes more difficult when the dimension increases, as
These phenomena are respectively denoted as masking and swamping, we will cover it in details during the robust statistics module of this course!
Let us conclude this lab session by looking at some DD-plots. For two probability distributions \(F\) and \(G\) , both in \(R ^d\), and \(D(\cdot)\) an affine-invariant depth, we can define depth vs. depth plot being very useful generalization of the one dimensional quantile-quantile plot:
\[ D D(F, G)=\left\{\left(D_{F}(x), D_{G}(x)\right) \text { for all } x \in \mathbb{R}^{d}\right\} \]
Its sample counterpart calculated for two samples \(\mathbf{X} =\left\{X_1, \cdots, X_n \right\}\) from \(F\), and \(\mathbf{Y} =\left\{Y_1, \cdots, Y_m \right\}\) from \(G\) is defined as
\[ D D\left(F_{n}, G_{m}\right)=\left\{\left(D_{F_{n}}(z), D_{G_{m}}(z)\right), z \in\{\mathbf{X} \cup \mathbf{Y}\}\right\} \]
The ddPlot function automatically computes and plots
it:
df_good <- df_3[1:n_good,]
df_out <- df_3[(n_good+1):n,]
ddPlot(x = df_good,y = df_out,depth_params = list(method='Tukey'))
## DDPlot
##
## Depth Metohod:
## Tukey
They would line up if they came from the same distribution
It is easy to manually build a DD-plot
depth_good <- depth(u = df_3,X = df_good,method = "Tukey")
depth_out <- depth(u = df_3,X = df_out,method = "Tukey")
plot(depth_good,depth_out, col="blue", xlab="X depth", ylab="Y depth", main= "Depth vs. depth plot")
grid(10, 10, col="grey50", lty=1)
abline(0,1, col="grey50")
Clearly the distributions are not the same… What if we had some extra samples coming from \(F\) (e.g., the trivariate normal centered at \(0\))?
n_extra <- 100
df_extra <- data.frame(mvrnorm(n_extra, mu_good, sigma_common))
ddPlot(x = df_extra, df_good,depth_params=list(method='Tukey'))
## DDPlot
##
## Depth Metohod:
## Tukey
Indeed this time we can conclude the same distribution generated the data.
Breakdown point
fracs <- seq(0.1, .5, length.out=25)
set.seed(41703192)
n_good <- 100
df.clean <- data.frame(mvrnorm(n_good, mu_good, sigma_common))
multivar.mean.clean <- colMeans(df.clean)
multivar.median.clean <- depthMedian(df.clean,depth_params = list(method='Tukey'))
mean.diffs <- numeric(length(fracs))
median.diffs <- numeric(length(fracs))
for (i in 1:length(fracs)){
n=100
# sample points
frac <- fracs[i]
n_good=ceiling(n*(1-frac))
n_outliers=n-n_good
df_contaminated_normals = data.frame(rbind(
as.matrix(df.clean[1:n_good, ]),
mvrnorm(n_outliers, mu_outliers, sigma_common)
))
current.mean <- colMeans(df_contaminated_normals)
current.median <- depthMedian(df_contaminated_normals, depth_params = list(method='Tukey'))
mean.diffs[i] <- max(abs(current.mean - multivar.mean.clean) )
median.diffs[i] <- max(abs(current.median - multivar.median.clean))
}
plot(x=fracs, y=mean.diffs, main="Empirical Breakdown point")
lines(x=fracs, y=median.diffs, col="red")
sim.data <- rnorm(20)
stat1 <- function(data){base::mean(data)}
stat2 <- function(data){median(data)}
L = 25
grid.for.x <- seq(-1.5*min(sim.data), max(sim.data)*1.5, length.out=L)
stat1.0 <- base::mean(sim.data)
stat2.0 <- stat2(sim.data)
EIF1 <- numeric(L)
EIF2 <- numeric(L)
for (i in 1:length(grid.for.x)){
sim.data.i <- sim.data
sim.data.i[1] <- grid.for.x[i]
EIF1[i] <- abs(base::mean(sim.data.i) - stat1.0)
EIF2[i] <- abs(stat2(sim.data.i) - stat2.0)
}
plot(x=grid.for.x, y=EIF1, ylim=c(-.05, .12))
lines(x=grid.for.x, y=EIF2, col="red")
An alternative is the empirical influence value: the
(estimated) the derivative of a sample estimator with respect to each
statistical unit in the sample of its input. One must be careful on the
definition of the statistic. It is implemented in the boot
package and the statistic (i.e. a functional of the sample), has to be
an (R) function of the data matrix and a vector of
weights, i.e the weight associated to each statistical unit
(row in the matrix), see more information here
For simplicity, we try it with the \(t\) statistic on univariate data: \(t(\mathbb{X}) = \frac{\hat{\mu}(\mathbb{X})-\mu_0}{\hat{\sigma}(\mathbb{X})}\) where \(\mathbb{X}\) denotes the data sample, with \(hat{\mu}(\mathbb{X})\) denoting the sample mean, \(\mu_0\) a hypothesised mean and \(\hat{\sigma}(\mathbb{X})\) the sample standard deviation
library(boot)
##
## Caricamento pacchetto: 'boot'
## Il seguente oggetto è mascherato da 'package:robustbase':
##
## salinity
mean.statistic <- function(data, w){
return ()
}